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Abstract 



We study the thermodynamics and the cohapse of a seh-gravitating gas of Brown- 
ian particles. We consider a high friction Umit in order to simpUfy the problem. This 
C^ . results in the Smoluchowski-Poisson system. Below a critical energy or below a critical 

^ I temperature, there is no equilibrium state and the system develops a self-similar collapse 

— Jj ' leading to a finite time singularity. In the micro canonical ensemble, this corresponds to 

P^ ■ a "gravothermal catastrophe" and in the canonical ensemble to an "isothermal collapse" . 

O I Self-similar solutions are investigated analytically and numerically. 

^ : 1 Introduction 
in 

^ . The thermodynamics of self-gravitating systems displays intriguing features due to the exis- 
f^^ ! fence of negative specific heats, inequivalence of statistical ensembles and phase transitions 
O I associated with gravitational collapse [|l|. Thermodynamical equilibrium of a self-gravitating 
system enclosed within a box exists only above a critical energy Ec = —0.335GM^/R or above 
a critical temperature T^. = GMm/2.52kR and is at most a metastable state, i.e. a local max- 
imum of a relevant thermodynamical potential (the entropy in the microcanonical ensemble 
and the free energy in the canonical ensemble) [^ ^. For T < T^ or E < E^, the system 
'^ I is expected to collapse. This is called "gravothermal catastrophe" or "Antonov instability" 
Q I in the microcanonical ensemble (MCE) and "isothermal collapse" in the canonical ensemble 
O ■ (CE). Dynamical models appropriate to star formation 0] or globular clusters f^, ||, 0, H show 
that the collapse is self-similar and leads to a finite time singularity (i.e., the central density 
^ ! becomes infinite in a finite time). The value of the scaling exponent in the density profile 
p ~ r~" depends whether the system evolves at fixed temperature (in which case a = 2 results 
from dimensional analysis) or if its temperature is free to diverge (in which case the value of 
the exponent is non trivial and often close to 2.2). It is found in general that the shrinking of 
the core is so rapid that the core mass goes to zero at the collapse time although the central 
density is infinite. 

In this paper, we introduce a simple model of gravitational dynamics which exhibits similar 
features and which can be studied in great detail. Specifically, we consider a gas of self- 
gravitating Brownian particles enclosed within a spherical box. For simplicity, we take a high 
friction limit and reduce the problem to the study of the Smoluchowski-Poisson system. In the 
simplest formulation, the temperature is constant (canonical description). We also consider 



the case of an isolated medium with an infinitely large thermal conductivity so that its tem- 
perature is uniform in space but varies with time in order to conserve energy (microcanonical 
description). The interest of these models is their relative simplicity which allows for a complete 
theoretical analysis, while keeping all the richness of the thermodynamical problem: inequiv- 
alence of statistical ensembles, phase transitions, gravitational collapse, finite time singularity, 
persistence of metastable states, basin of attraction... These models are consistent with the 
first and second principles of thermodynamics and give a dynamical picture of what happens 
when no equilibrium state exists. However, in view of their considerable simplification, it is 
not clear whether these models can have astrophysical applications although connections with 
the dynamics of dust particles in the solar nebula and the process of "violent relaxation" in 
coUisionless stellar systems are mentioned. 

The paper is organized as follows. In Sec. |2.1| , we introduce the Smoluchowski-Poisson 
(SP) system for a gas of self-gravitating Brownian particles and list its main properties. In 
particular, we make contact with thermodynamics and show that the SP system satisfies a 
form of iiT-theorem. In Sec. p.2| , we discuss the existence of stationary solutions of the SP 



system and the relation with maximum entropy states. In Sec. p.3| , we perform a linear 
stability analysis of the SP system. We show that a stationary solution is linearly stable if and 
only if it is a local entropy maximum and that the eigenvalue problem for linear stability is 
connected to the eigenvalue problem for the second order variations of entropy studied in Refs. 
[^, [T^. In Sec. 1^, we consider the case of gravitational collapse and exhibit self-similar solutions 



of the SP system. Since the particles are confined within a box, there is a small deviation to 
the purely self-similar regime and we describe this correction in detail. 

In Sec. H, we perform various numerical simulations of the SP system for different initial 
conditions. We check the results of thermodynamics, namely the existence of equilibrium states 
for E > Ec and T > T^. and the gravitational collapse otherwise. We find that the collapse 
proceeds self-similarly with explosion, in a finite time tcoiu of the central density while the 
core radius shrinks to zero. In MCE this is accompanied by a divergence of temperature and 
entropy. In the limit t — > tcoiu "we find the scaling laws pqTq ~ 1 and p/po ~ ('"/'"o)~"- The 
scaling exponent is a = 2 in CE and a ~ 2.21 in MCE. In CE, the invariant profile p/po = 
/{r/ro) can be determined analytically. The collapse time diverges like tcou ~ {Ec — E)"^^"^ 
and tcoii ~ {Tc — T)^^/^ as we approach the critical energy Ec and critical temperature T^. We 
also study the linear development of the instability (for unstable isothermal spheres) and show 
that the density perturbation Sp/p presents several oscillations depending on the value of the 
density contrast. In particular, at the points of marginal stability in the series of equilibria, the 
perturbation Sp/p has a "core-halo" structure in the microcanonical ensemble but not in the 
canonical ensemble in agreement with theory [^, |T^. 

2 Self-gravitating Brownian particles 

2.1 The Smoluchowski-Poisson system 

We consider a system of small particles with mass m immersed in a fiuid. We assume that the 
fluid imposes to the particles a friction force — ^v and a stochastic force R(t). This random force 
may mimic ordinary Brownian motion (i.e. the collisions of the fiuid particles onto the solid 
particles) or fiuid turbulence. We assume in addition that the particles interact gravitationally 
with each other. Therefore, the stochastic Langevin equation describing the motion of a particle 



reads 

(1) ^ = -ev + F(r,t) + R(t), 

where F = — V^ is the gravitational force acting on the particle. For simplicity, we shall 
assume that the stochastic force is delta-correlated in time and set 

(2) {R{t)R{t')) = 6D 6{t-t'), 

where D measures the noise strength of the Langevin force. In order to recover the Maxwell- 
Boltzmann distribution 

at equilibrium, the diffusion coefficient and the friction term must be related according to the 



Einstein relation D = C,T. Applying standard methods [11|, we can immediately write down 



the Fokker-Planck equation associated with this stochastic process: 

dt dr dv dv\ \dv J ) 

This is the familiar Kramers equation but, when self-gravity is taken into account, it must be 
coupled to the Poisson equation 

(5) A$ = AttGp, 

where G is the gravitational constant. This makes its study much more complicated than usual. 
The Kramers-Poisson (KP) system was first introduced in astrophysics by Chandrasekhar [O 
in his stochastic theory of stellar dynamics (for, e.g., globular clusters). In that context, the dif- 
fusion and the friction arise self-consistently as the result of the fiuctuations of the gravitational 
field. An equation of the form (H) was also proposed as an effective dynamics of collisionless 



stellar systems (on a coarse-grained scale) during the period of violent relaxation ||T3|, |14 . 

In order to simplify the problem, in a first approach, we consider a high friction limit 
^ —^ +00. Then, it is possible to neglect the inertial term in the Langevin equation (0). The 
Fokker-Planck equation describing this high friction limit is the Smoluchowski equation 

(6) |^ = v|i(TVp + pV<l>) 

with a diffusion coefficient D' = T/C, and a drift term proportional to the gravitational force. 
The ordinary Smoluchowski equation describes the sedimentation of colloidal suspensions in an 
external gravitational field. Since it is a prototype of kinetic equations, it is clearly of great 
interest to consider the extension of this model to the case where the potential is not fixed but 
related to the density of the particles via a Poisson equation, like in the gravitational case. 

The Smoluchowski equation can be interpreted equivalently as a continuity equation for the 
density p with a velocity field 

(7) u = -V-Vp + V<l> 

where —TVp is the pressure force and — pV$ the gravitational force. At equilibrium, the two 
terms balance each other and the Boltzmann distribution (|^) establishes itself. Physically, the 



high friction hmit supposes that there are two time scales in the problem. On a short time scale 
of the order of the friction time ^~^ ^ tdyn, the system thermalizes and the distribution function 
becomes Maxwellian with temperature T (this is obvious if we take the limit D = ^T -^ +00 in 
the r.h.s. of Eq. (Q)). Then, on a longer time scale of the order of the dynamical time t^yn, the 
particle distribution p(r, t) tends to evolve towards a state of mechanical equilibrium described 
by the Boltzmann distribution (^). Note that the opposite assumptions are made for globular 
clusters [§|, |^, 0: the system is assumed to be in mechanical equilibrium and the evolution is 
due to thermal transfers between the core and the halo. Our model of self-gravitating Brownian 
particles could find applications for the dynamics of dust particles in the solar nebula and the 
formation of planetesimals by gravitational instability (see, e.g., Ref. [|15|). In that context. 



the dust particles experience a friction with the gas modeled by Stokes or Epstein's laws and 
the high friction limit may be relevant. On the other hand, the diffusion of the particles could 
result from a stochastic component of the force or from fluid turbulence. This would be just a 
first approach because the physics of planetesimal formation is more involved than our simple 
model. 

Since the system described previously is in contact with a heat bath, the proper statistical 
treatment is the canonical ensemble in which the temperature T is fixed. In order to test 
dynamically the inequivalence of statistical ensembles for self-gravitating systems, we would 
like to introduce a simple model corresponding to the microcanonical ensemble, i.e. with strict 
conservation of energy E. In fact, when a Brownian particle moves with its terminal velocity in 
a gravitational field, the work of the force ought to be converted into heat. If the medium acts 
as a thermostat with an infinite volume and with rapid dissipation of heat, we can disregard 
the variation of temperature and we get the isothermal model discussed previously. However, 
if we are to keep track of local heating, the temperature will depend on space and time and we 
need to set up a model in which energy is conserved. Such a generalization of Brownian theory 



has recently been developed by Streater |T^ in the case of an external gravitational potential. 
This statistical dynamics approach [|r^ leads to coupled nonlinear equations for the density 
p{r,t) and the temperature T{r,t) which are consistent with the first and second principles of 
thermodynamics. Such equations can be derived from a microscopic model involving Brownian 
particles and heat particles modeled as quantum oscillators. A generalization of these equations 
for self-gravitating Brownian particles has been proposed by Biler et al. [^ . It consists of the 
Smoluchowski-Poisson system (^ (|^) coupled to a diffusion equation for the temperature 

(8) --(pT) = V(AVT) - V(TJ) - JV$, 

where J is the diffusion current in Eq. (P). However, this model still remains complicated for 
a first approach. Since our main purpose is to illustrate in the simplest way the basic features 
of the thermodynamics of self-gravitating systems (inequivalence of ensembles, gravothermal 
catastrophe, isothermal collapse, phase transitions, basin of attraction...), we shall consider an 
additional approximation and let the thermal conductivity A in Eq. (H) go to +cx). In that 
case, the temperature is uniform but still evolving with time according to the law of energy 
conservation (first principle): 

(9) E=^-MT{t) + ]^j p^dh. 

The first term in the r.h.s is the kinetic energy K = J f'^d^rd^^ for a Maxwellian distribution 
function with temperature T (local thermodynamical equilibrium) and the second term is the 
gravitational energy of interaction. Equations (H) (Bl) (pi) lead to a simple microcanonical 



model for self-gravitating systems with a lot of attractive properties. The Cauchy problem for 



this system of equations was studied by Rosier [T^]. These equations were first proposed by 



Chavanis et al. [^ as a simplified model of "violent relaxation" by which a stellar system 



initially far from mechanical equilibrium tries to reach an isothermal state on a few dynamical 



times [jT3|, ^. In that context, the engine of the evolution is the competition between pressure 
and gravity, like in Eq. (^. This particular equation corresponds to an overdamped evolution 
but more general equations taking into account inertial terms are also proposed in Ref. |^ . 

It is easy to show that the SP system admits a form of if -theorem for an appropriate 
thermodynamical potential (second principle). The microcanonical ensemble is characterized 
by the specification of mass M and energy E. The thermodynamical potential is the entropy 

(10) S= -M+-M\n{2nT)- I plnpdh, 

which is the form of the classical Boltzmann entropy S = — j fin fd^rd^w for a Maxwellian 
distribution function with temperature T. Then, it is easy to show, using Eqs. (H) and (||) that 
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(11) S = j ^{TVp + pV*)^ dh > 0. 

Therefore, the entropy plays the role of a Lyapunov function for our microcanonical model. 
The canonical ensemble is characterized by the specification of mass M and temperature T. It 
is straightforward to show that the SP system (|^) satisfies a relation similar to Eq. ( |TID for 
the free energy (more precisely the Massieu function) J = S — [3E. It can be noted that the 
Kramers equation @) and the Smoluchowski equation (|^) can also be derived from a variational 
formulation [^, called the Maximum Entropy Production Principle (M.E.P.P.). This makes 
a direct relation between the dynamics and the thermodynamics. Since the SP system with 
the constraint (^ obeys the same conservation laws and H-theorem as more realistic models 
such as Landau- Poisson system and coarse-grained Vlasov-Poisson system |1^, it should 



exhibit qualitatively similar properties even if the details of the evolution are expected to differ 
in many respects. 

To properly define our system of equations, we must specify the boundary conditions. We 
shall assume that the system is non rotating and restrict ourselves to spherically symmetric 
solutions. In addition, we shall work in a spherical box of radius R to avoid the well-known 
infinite mass problem associated with isothermal configurations. In that case, the boundary 
conditions are: 

(12) ^(0) = 0. *W = -^, T^ + P^-0. 

The first condition expresses the fact that the gravitational force at the center of a spherically 
symmetric system is zero. The second condition defines the gauge constant in the gravitational 
potential. Finally, the last condition insures that the total mass is conserved (we have used the 
Gauss theorem dr^ = GM/r'^ to simplify its expression). 

For spherically symmetric systems, it is possible to reduce the SP system to a single partial 
differential equation for the mass profile M{r, t) = An J^ pr "^dr' . Multiplying both sides of Eq. 
(|^) by r^ and integrating from to r we obtain after straightforward algebra 

DM s_ir^5^M 2TdM GM{r,t)dM 

dt ^ ^\ dr"^ ' r dr '' ^2 g^ ^ 



The appropriate boundary conditions are now M(0, t) = and M{R, t) = M. The potential 
energy can be expressed in terms of M{r, t) as 



(14) W = - \^ {r,t)dr. 

Jo r dr 

It is possible to simplify Eq. (ITSp a little more by introducing the new coordinate m = r^ so 



that 

(15) ^^^u,t) = 9Tu'/'^{u,t) + 3GM{u,t)^{u,t). 

Finally, we note that the Krammers-Poisson (KP) system satisfies a form of Virial theorem: 

(16) l«+l^|.2A- + W-3p.V. 

where I = J pr'^d^r is the moment of inertia (we have properly taken into account the pressure 
on the box). The difference with the usual Virial theorem is the occurrence of a damping term 
1^/ due to friction. In the high friction limit, we get 

(IT) l|=l(2A- + ,._3p,l.). 

This expression can also be directly obtained from the SP system. 

2.2 Stationary solutions and maximum entropy states 

The stationary solutions of the SP system are given by the Boltzmann distribution (|^) in which 
the gravitational potential appears explicitly. The Boltzmann distribution can also be obtained 
by maximizing the entropy S at fixed mass and energy or by maximizing the free energy J at 
fixed mass and temperature. The gravitational potential is determined self-consistently by 
solving the mean field equation 

(18) A$ = AirGAe-^"^, 

obtained by substituting the density (^ in the Poisson equation (^. This Boltzmann- Poisson 
equation has been studied in relation with the structure of isothermal stellar cores P3[ and 



globular clusters p2| . It is well-known that the density of an isothermal gas decreases at large 
distances like r~^ resulting in the infinite mass problem if the system is not bounded. 

The equilibrium phase diagram {E, T) of isothermal configurations confined within a box is 
represented in Fig. [l| where we have plotted the normalized inverse temperature rj = (3GM/R as 
a function of the normalized energy A = —ER/GM'^. The curve has a striking spiral behavior 
parameterized by the density contrast TZ = p{0)/p{R) going from 1 (homogeneous system) 
to +00 (singular sphere) as we proceed along the spiral. There is no equilibrium state above 
Ac = 0.335 or rjc = 2.52. In that case, the system is expected to collapse indefinitely. It is also 
important to recall that the statistical ensembles are not interchangeable for systems with long- 
range interaction, like gravity. In the microcanonical ensemble, the series of equilibria becomes 
unstable after the first turning point of energy (MCE) corresponding to a density contrast of 
709. At that point, the isothermal spheres pass from local entropy maxima to saddle points. 
In the canonical ensemble, the series of equilibria becomes unstable after the first turning 
point of temperature (CE) corresponding to a density contrast of 32.1. At that point, the 
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Figure 1: Equilibrium phase diagram for classical isothermal spheres. The spiral rolls up 
indefinitely towards the singular isothermal sphere ps = l/2TiG(3r'^. 



isothermal spheres pass from maxima of free energy to saddle points. It can be noticed that 
the region of negative specific heats between {CE) and {MCE) is stable in the microcanonical 
ensemble but unstable in the canonical ensemble as expected on general physical grounds |l|. 
The thermodynamical stability of isothermal spheres can be deduced from the topology of the 
j3 — E curve by using the method of Katz [^ who has extended Poincare's theory of linear 
series of equilibria. The stability problem can also be reduced to the study of an eigenvalue 
equation associated with the second order variations of entropy or free energy as studied by 
Padmanabhan [§ in MCE and Chavanis [|l^ in CE. The same stability limits as Katz are 
obtained but this method provides in addition the form of the density perturbation profiles 
that trigger the instability at the critical points. We also recall that isothermal spheres are at 
most metastable: there is no global maximum of entropy or free energy for a classical system 
of point masses in gravitational interaction p[|. 



2.3 Linear stability analysis 

We now perform a linear stability analysis of the SP system. Let p, T and $ refer to a stationary 
solution of Eq. (|^) and consider a small perturbation 5p, 5T and ^<l> around this solution that 
does not change energy and mass. Since a stationary solution of the SP system is a critical point 
of entropy, we must assume A < Ac for a solution to exist. Writing 5p ~ e'^* and expanding 
Eq. (^ to first order, we find that 

(19) 






dp ^dSp (i$ (i(5$ 



^ \ dr dr 



+ ^p-r^ P' A 

dr dr 



It is convenient to introduce the notation 



(20) 



5p 



1 dq 



Anr^ dr 

Physically, q represents the mass perturbation q{r) = 5M{r) = J^ Anr' 5p{r')dr' within the 
sphere of radius r. It satisfies therefore the boundary conditions g(0) = q{R) = 0. Substituting 



Eq. (^) in Eq. ([T9|) and integrating, we obtain 

(21) ^g = 4nST^P + T^(Y4)+V4f + inf.'"'' 



r^ dr dr\r'^drj r'^ dr dr dr 

where we have used g(0) = to ehminate the constant of integration. Using the condition of 
hydrostatic equihbrium Tdp/dr + pd^/dr = and the Gauss theorem d6^/dr = Gq/r"^, we can 
rewrite Eq. (^) as 

A^ 6Td(l> 1 d f 1 dq\ 1 Idqdp Gq 

(22) 7— 7^3? = -7^-z: + 7ir:3i z?x: -7r3z?x:x: + 



AnpTr^ T^ dr Anp dr \r'^ dr J Airp"^ r^ dr dr Tr"^ ' 

or, alternatively, 

d ( 1 dq\ Gq \i bTd^ 



dryAnpr'^drJ Tr"^ AixpTr'^ T"^ dr 

From the energy constraint (^ we find that 

2 1"^ 2 f^ d^ 
(24) 6T = / 5p^Anr^dr = / q— dr. 

Hence, our linear stability analysis leads to the eigenvalue equation 

(25) 

where 

(26) 

where we recall that g(0) = q{R) = 0. Eq. (^) is similar to the eigenvalue equation associated 
with the second order variations of entropy found by Padmanabhan p. In particular, they 
coincide for marginal stability (A = 0). More generally, it is proven in Appendix that a 
stationary solution of Eq. (|^) is linearly stable if and only if it is a local entropy maximum. 
The zero eigenvalue equation was solved by Padmanabhan ||^. It is found that marginal stability 
occurs at the point of minimum energy A = Ac, in agreement with Katz p4| approach, and that 



d f 1 dq\ Gq 2V d$ 
dr\47rpr^dr) ' Tr^ SMT'^ dr 


47rpTr2 







the perturbation 6p/p that induces instability (technically the eigenfunction associated with 
A = 0) has a "core-halo" structure (i.e., two nodes). It is also argued qualitatively that the 
number of oscillations in the profile Sp/p increases as we proceed along the series of equilibria, 
see Fig. |I], up to the singular sphere (i.e for higher and higher density contrasts). Of course, on 
the upper branch of Fig. |I|, the eigenvalues A are all negative (meaning stability) while more 
and more eigenvalues become positive (meaning instability) as we spiral inward for TZ > 709. 

If we fix the temperature T instead of the energy E, the eigenvalue equation becomes (take 
6T = 0m Eq. (ESl)): 



('27) ^ ( ^ ^^^ + '^'^ ^^ 



dr yAiT pr'^ dr J Tr"^ AirpTr'^ 



This is similar to the equation obtained by Chavanis [|10| by analyzing the second order varia- 
tions of free energy. The case of marginal stability (A = 0) coincides with the point of minimum 



temperature rj = rjc like in Katz |23] analysis. It is found that the perturbation 6p/p that in- 
duces instability at rj = rjc in the canonical ensemble has not a "core-halo" structure (it has 
only one node). 
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3 Self-similar solutions of the Smoluchowski-Poisson sys- 
tem 

3.1 Formulation of the general problem 

We now describe the collapse regime and look for self-similar solutions of the SP system. 
Restricting ourselves to spherically symmetric solutions and using the Gauss theorem, we obtain 
the integrodifferential equation 

We look for self-similar solutions in the form 

r \ I 1 ^ 



(29) ,(.,() =p„(()/^_J, ,, ^^^^ 

where the density po(^) is of the same order as the central density p(0, t) and the radius tq is of 
the same order as the King radius r^ = (9T/47rGp(0))^/^ which gives a good estimate of the 



core radius of a stellar system p2[. Substituting the ansatz (p9D into Eq. (p8|), we find that 



(») >^^ - ^lir^^'f^' - f^M^i^''^' ^ i/(.)f /(.')ww)}. 

where we have set x = t/tq. The variables of position and time separate provided that there 
exists a such that pqTq ~ 1. In that case, Eq. ( ^0]) reduces to 



Assuming that such a scaling exists implies that {^/GpQ){dpo/dt) is a constant that we arbi- 
trarily set to be equal to 1. This leads to 

(32) Po{t) = ^it,M-t)-\ 

so that the central density becomes infinite in a finite time tcou while the core shrinks to zero as 
'^o ~ {tcoii — ty^"- Since the collapse time appears as an integration constant, its precise value 
cannot be explicitly determined. The scaling equation now reads 

(33) fix) + -xfix) = \^L'(f'{x) + l/(x) r f{x')47rx''dx' 

a x'' ax { \ s^ Jq 

which determines the invariant profile f{x). Alternative forms of Eq. (|33|) are given in Appendix 
^ If one knows the value of a, Eq. ( |5BD leads to a "shooting problem" where the value of /(O) 
is uniquely selected by the requirement of a reasonable behavior for f{x) at large distances (see 
below). As f{x) — i> for large x, we can only keep the leading terms in Eq. (|33D , which leads 
to f{x) ~ x^°' when x -^ +oo. 

The velocity profile defined by Eq. (|^) can be written 

(34) u{r,t) = -voit)V' 



ro{t)J' 



with 



(35) vo{t) = ^ and V{x) = Q^ + ^ T f{x')x' dx' . 

The invariant profile V{x) has the asymptotic behaviors V{x) ~ x when x —^ and V{x) ~ 1/x 
when X — >■ +cxd. On the other hand, the mass profile can be written 

(36) M{r,t) = Mo{t)gi '' 



with 

(37) Mo(t) = po^o and g{x) = An f{x')x'^dx' . 

Jo 

The invariant profile g{x) has the asymptotic behaviors g{x) ~ x^ when a; — > and g{x) 
when X — » +oo. 



~a;3-" 



3.2 Canonical ensemble 

In the canonical ensemble in which the temperature T is a constant, Eq. (pOj) leads to a = 2 
(the particular case T = is treated in Appendix [B]). In that case, the scaling equation (|33| ) 
can be solved analytically (see Appendix ^) and the invariant profile is exactly given by 

(38) fix) = -LJ±4^. 

47r (1 + ^)2 

This solution satisfies /(O) = ^ and f(x) ~ ^ as x — > +cxd. From Eq. (|32D, we find that the 
central density and the core radius evolve with time as 

(39) p(0, t) = po(t)/(0) = ^(tcou - tr\ ro(t) = (|) (t,„, - t)^/^ 



On the other hand, using Eq. (|38|), we find that the velocity profile and the mass profile are 



given by Eqs. (|3^ ) and ( |36D with 

(40) ^;o(t)=f^) ^t,oii-t)-'/^ and ^(x) 



-^\ / \l/0 -. -,^/\ '^^ 



6 + x2 



T^V^\ .,/. , _ 4x3 



(41) Mo(t)= -— (tcoii-t)^/' and 



^ia;j 



^^GV ^ '"" ^ "^ ^ 2 + x2 

At t = tcoii, the scaling solutions (^) (|40|) and (^TJ) converge to the singular profiles 

T 2T 4T 

(42) p(r, t = tcoii) = —pT^, u{r, t = t^ou) = -7—, M{r, t = t^oii) = -pr^- 

vrGr^ 4r G 

It is interesting to note that the density profile (^2]) has the same r-dependence as that of the 
singular solution to the static isothermal gas sphere p = l/27rG/5r^ [^, the two profiles just 
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differing by a factor of 2. Tlierefore, tlie relationsliip between tlie density and tlie gravitational 
potential in the tail of the scaling profile is given by a Boltzmann distribution 



(43) p ~ Ae 



2T 



with a temperature 2T instead of T. A r^^ decay of the density at large distances was also 
found by Penston Q] in his investigation of the self-similar collapse of isothermal gas spheres 
described by the Euler equations. This is a general characteristic of the collapse in the canonical 
ensemble (T = Cst.). It should be noticed that the free energy does not diverge at tcou although 
the system undergoes a complete collapse. Therefore, at t = tcoii, the density profile is not a 
Dirac peak contrary to what might be expected from rigorous results of statistical mechanics 



25|. In fact, there is no contradiction because the Dirac peak is formed during the post collapse 



evolution 

We now show that the self-similar solution (^9]) is not sufficient to quantitatively describe 
the full density profile (especially when r ~ i?). To understand the problem, let us calculate 
the mass contained in the scaling profile at t = tcoii- Using Eq. (|42|), we have 

f^ T . ^ AR 

(44) M...„, = y^ ^^4^^^^ = G^- 

The mass M scaling is finite but, in general, it is not equal to the total mass M imposed by 
the initial condition. This means that there must be a non-scaling contribution to the density 
which should contain the remaining mass (possibly negative when M < Mgcaiing, i-e. r/ < 4). 
That the scaling solution ( p9D is not an exact solution of our problem is also visible from the 
boundary conditions. Indeed, according to Eq. ([l^ ) we should have 

,^ , d\np [3GM 

This relation is clearly not satisfied by Eq. (^) except for the particular value t] = 2. These 
problems originate because we work in a finite container. The scaling solution ( PU[ ) would 
be exact in an infinite domain but, in that case, the total mass of the system is infinite. In 
addition, if we remove the box, the isothermal spheres are always unstable and the interesting 
bifurcations between equilibrium and collapsing states are lost. 

Strictly speaking, we expect that the self-similar solution (|29|) will describe the density 
behavior in the scaling limit defined by 

(46) t — ^ tcoii or ro — ^ 0, and x = r /tq fixed. 

For the reasons indicated above, it probably does not reproduce the density near the edge of 
the box, that is for r ~ i? ^ tq. Therefore, we write another equation for the density, making 
the following ansatz: 

(47) ,(.,,) = ,„(,)/(_Z_)+_I_^(,,,), 

where F{r,t) is the profile that contains the excess or deficit of mass. For t = tcoii, we have 

(48) ^^"'^-"^ = 4^(^ + ^^"0' 

and it would be desirable to find an approximate expression for the function F{r) = F{r,tcoii)- 
A differential equation for F{r) can be obtained by substituting the ansatz (^) in the dynamical 
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equation (pSD and taking the limit t = tcou- We need first to discuss the term dp/dt{r,tcoii)- 
For t -^ tcoii, we can use the expansion of the function /(x), given by Eq. (^), to second order 
in 1/x^ to get 



(49) 



p{r,t) 



Poro 






T 

4^ 



F(r,t). 



Then, using Eqs. ( p9|) and (|3^) , we obtain to first order in tcoii — t: 

T2 



m 

leading to 
(51) 



p{r,t) = p{r,tcou) + 






AttG^ 



T2 



^4 rp g^ {^^ ^coll) 



(tcoU — t) + ... 



47rG^ 



__8_ ^d£ 
r^~^ T dt '^'''^™"'' 



The trouble is that we do not know the function dF/dt{r,tcoii)- It is possible, however, to 
derive an exact integral equation that it must satisfy. Since the exact profile p{r, t) conserves 
mass, we have just before tcoif- 



(52) 



R 






^(rA,,)r'^dr = 0. 



The scaling profile Pscaiing{i^,t) is an exact solution of Eq. (p8D but it does not conserve mass. 
Multiplying Eq. (^) by r^ and integrating from r = to i?, we get 



(53) / 
Jo 



^ P scaling 



dt 



r,-t,Mydr 



-K I rp^ P scaling t -p\ '^ ^''J- scaling \ 

~Y\ 'dr P scaling ^^ I 



1^ ''coli 



2T^ 
TiGiR' 



where we have used Eqs. ( ^2)) and (|4^) to obtain the last equality. Now, subtracting Eqs. (|5^ ) 
and (|53|), using Eq. ( ^Tj ) and passing to the limit t -^ tcoii, we find that 



(54) 



R 



dF , , 2 , 8T 

-T^{f,tcoii)r dr = ——;-. 

dt ^ ' ^ ^i? 



This relation implies in particular that we cannot take {dF/dt){r, tcou] 
it is likely that F{r,t) involves combinations of the type 



OinEq. (|l|). In fact. 



(55) 



F{r,t)r^Pof{r/royF{r), l^(r^ + crl)F{r), F{^r^ + cr^), ... 



which reduce to F{r) in the limit t —>■ tcoii- Considering the time derivative of these expressions 
at t = tcoii, we find that they take only one of the two forms F{r)/r'^ and F'{r)/r. We are 
therefore led to make the following ansatz : 



(56) 






Fir) F'ir) 

fy-^ ly 



where a and b are some unknown constants which will be determined by an optimization 
procedure (see below). If we substitute the ansatz (|47|) in Eq. (pSf), take the limit t = tcoii and 
use Eqs. (0) and (|56D , we find after some simplifications that F{r) satisfies the differential 
equation 



(57) 



r^F" + (6 - b)rF' + r^F^ + (8 - a)F + F' [ 

Jo 



F(x)x'^dx 



/ F{x)x^dx = 0. 
Jo 
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Interestingly, the final profile equation (^) is not obtained by setting dp/dt = in the dynam- 
ical equation as, even in the stationary looking tail, dp/dt is in fact of order 1 due to the fast 
collapse dynamics. 

Equation (^1\ ) leads to another "shooting problem", starting this time from r = R. The 
value F{R) is selected by imposing the condition that the total mass is M. This yields 



"R / /I D\ 

where AR/PG = M scaling is the mass included in the scaling part. Moreover, F\R) is fully 



(58) / F{rydr = (3g(M-^\ 

Jo \ P^y 



determined by the boundary condition (12) at r = i? which implies, together with Eq. 



Finally, the exact relation ( [5^ ) combined with Eq. ( pBD imposes the condition 

F(r)dr^hRF(m = - . 

R 



(60) (a -6) / F{r)dr + bRF{R) 

Jo 



In order to determine the values of a and b we shall require that the value of the total density 
at r = i? is maximum, as the system would certainly tend to expel some mass if it were not 
bound to a sphere (recall that the profile F arises because of boundary effects). In addition, 
Eq. (|60|) implies that F is integrable, so that the optimization process should be performed 
including this constraint (if F is integrable, then Eq. ( |60D is automatically satisfied as it is 
equivalent to the conservation of mass). In the section devoted to numerical simulations, we 
study F numerically and compare it with the numerical profiles obtained by solving the SP 
system. 

3.3 Microcanonical ensemble 

If the temperature is not fixed but determined by the energy constraint (^, then the exponent 
a is not known a priori. However, we have solved Eq. ( |5B| ) numerically for different values of 
a and found that there is a maximum value for a above which Eq. (pSf ) does not have any 
physical solution. This value amax = 2.20973304... is close to that found by Lynden-Bell & 
Eggleton (and, to some extent, by Cohn [^] and Larson |Q) in their investigation on the 
gravitational collapse of globular clusters. The common point between these models is that 
the temperature is free to diverge so the scaling exponent a cannot be determined from simple 
dimensional analysis. However, the agreement on the value of a is probably coincidental since 
our model differs from the others in many respects. 

In the present case, amax is just an upper bound on a not a unique eigenvalue determined 
by the scaling equations like in Ref. [0] for example. However, this maximum value leads to 
the fastest divergence of the entropy and the temperature so it is expected to be selected by 
the dynamics (recall that the SP system is consistent with a maximum entropy production 
principle pOj). Indeed, the temperature and the entropy respectively diverge like 



(61) T{t) ~ (t,„, - t)-^, S{t) ~ -^(^^ln(t,„, - t). 

Za 

Note that these divergences are quite weak as the exponent involved is small "'"""''^ = 0.0949133.. 
For a = Umax, the value of /(O) selected by the shooting problem defined by Eq. (|33D is 
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/(O) = 5.178.... Therefore, the central density evolves with time as 

(62) pi0,t)=5.m...^{t,ou-t)-\ 

The coefficient in front of {tcoii ^ t)^^ is approximately 10 times larger than for a = 2 (see 
Eq. (^)). The density proffie at t = tcoii is equal to 

(63) p(r,t = tcoii) = —, 

where K is a constant which is not determined by the scaling theory. Using Eq. (|63D and the 
Gauss theorem, we find that the relation between p and $ in the tail of the self-similar profile 
is that of a polytrope: 

(64) p~ ($- C5t.)^, 

with index n = a/{a — 2) ~ 10.53 for a = amax- 

We now address the divergence of the potential energy which should match that of the 
temperature (or kinetic energy) in order to ensure energy conservation. After an integration 
by parts, the potential energy can be written 

(65) w^-—-—jmfS.. 
Then, using the Gauss theorem, we obtain 

(66) w = -^-^l -U p{r')A7rr"'dr'\ dr. 

If we assume that all the potential energy is in the scaling proffie, we get a contradiction since 

(67) Wscaiing{t = Uu)---^ j ^ (f ^"^^^^ dA dr r^ - J T^-""^ dr, 

converges for a < 5/2. Since the temperature diverges with time for a = amax, the total 
energy cannot be conserved. This would suggest that a = 2 like in the canonical ensemble. 
We cannot rigorously exclude this possibility but a value of a close to amax — 2.21 is more 
consistent with the numerical simulations (see Sec. ^ and leads to a larger increase of entropy 
(in agreement with the MEPP). If this value is correct, the divergence of the gravitational 
energy should originate from the non scaling part of the profile which also accommodates for 
the mass conservation. In the following, a possible scenario allowing for the gravitational energy 
to diverge is presented. 

Let us assume that there exists two length scales ri and r2 satisfying r^ <^ ri <^ r2 <^ R 
with ro,ri,r2 —* for t —>■ tcoii such that the mass between ri and r2 is of order 1. The 
physical picture that we have in mind is that this mass will progress towards the center of the 
domain and form a dense nucleus with larger and larger potential energy. We assume that for 
ri < r < r2 the density behaves as 



(68) p{r,t)^^ 



7— a 

rv 



r7 
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so that this functional form matches with the scahng profile for r ~ ri. If we impose that the 
total mass between ri and r2 is of order 1, we get 

(69) / r'^dr ~ 1, i.e. r2 ~ ''^i'"^ , 

which shows that r2 ^ ri since a < 3. Now, the contribution to the potential energy of the 
density between ri and r2 which is assumed to be the dominant part is 

(70) Wr^- r^fr ^r'2 dA ' dr ~ -r^^^-^V^^^ ~ _^-(-7)/(3-7)^ 

where we have used Eq. (^) to get the last equivalent. Since the divergence of the potential 
energy must compensate that of the kinetic term we must have —W ~ |MT ~ Tq"" where we 
have used Eqs. (^9]) to get the last equivalent. This relation implies that tq and ri are related 
to each other by 

(71) ri ~ Tq 

Now, imposing ri ^ tq leads to 7 < 2. Therefore, any value of 7 < 2 leads to the correct 
divergence of W within this scenario. Note that Eq. (|68|) may arise from the next correction to 
scaling of the form 

(72) p{r,t) = Po/(r/ro) +pJ/i(r/ro) + ..., 

with /i(x) ~ x~'^ for large x and 7 < 1 for the first term to be dominant in the scaling regime. 
Matching the large x behavior of Eqs. ( pBD and ([7^), we obtain 

(73) pX ~ r7-^ 
which is equivalent to 

(74) n ~ r^ 



a — 'y 



Since 7 < 1, this implies that ri ^ tq, as expected. More precisely, comparing with Eq. (|7TD , 
we have 

(75) _^, + («-2)(3-,)_ 

a 

and we check that the condition 7 < 2 is equivalent to 7 < 1. 

3.4 Analogy with critical phenomena 

In this section, we determine the domain of validity of the scaling regime by using an analogy 
with the theory of critical phenomena. For simplicity, we work in the canonical ensemble but 
we expect to get similar results in the microcanonical ensemble. For 77 = 6^^ = f3GM/R close 
to rjc, we define 

/7RN 1^ - Vc\ \0c-6\ ^ ^ 
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For rj = rjc the central density p(0,t) goes to a finite constant poo when t — > +00. Writing 
Sp(t) = Poo — p(0, t) and using Eq. (^), which is quadratic in p, we argue that, for t] < rjc, 5 pit) 
satisfies an equation of the form 

/77N d5p 6p G 2 

where r plays the role of a correlation time which is expected to diverge for rj = t]c leading to a 
slow (algebraic) convergence of 6p towards at the critical temperature. Actually, for rj = rjc, 
Eq. (^) yields 



-1 



(78) Sp ~ t 

Now, if we stand slightly above the critical point (77 > rjc), we expect this behavior to hold up to 
a time of order tcoii for which the perturbation term proportional to (1/^)(T — Tc)Ap(0, t) ~ — e 
is of the same order as dp/dt ~ — 1/t^. This yields 

(79) t,oii--e-'/'^{v-Vc)-'^'. 

By analogy with critical phenomena, it is natural to expect that r has the same behavior for 
rj < r]c- 

(80) r~(r/,-r/)-i/2. 

Therefore, for rj < rjc and according to Eq. ([77|), 6p{t) tends exponentially rapidly to the 
equilibrium value 

(81) Poo - p(0, t = +00) = I r-i ~ [r], - r^fl\ 



This relation is consistent with the results obtained in the equilibrium study [|T0|, where the 
exact result 



p(o) ^ r 8 / 7]\y^^ 

poo [vc-2 V VcJi 

is derived close to the critical point. 

Another interesting question concerns the extent of the scaling regime which we expect to 
be valid for tcoii — t < St ^ e" . To compute z/, we integrate the dynamical equation in the 
regime where the perturbation —eAp dominates: 

(83) I . -.Ap. 
leading to 

(84) p(0,t)~/ Pexp{k^et)dkr^rQ^exp{rQ^et), 

where we have introduced an upper momentum cut-off of order Tq^ to prevent the integral from 
diverging. Indeed, the Laplacian of p should become positive for r :^ tq as A(r~^) = 2r~^ > 0. 
Thus, for e ^ 1, we expect that the density will first saturate to poo for a long time of order 
tcoii [see Eq. (^)], before rapidly increasing [see Eq. (|8l) ], and ultimately reaching the scaling 
regime [see Eq. (|29|)]. Comparing Eq. ( ^4]) with the density in the scaling regime p(0,t) ~ Tq'^, 
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we find that the scahng regime is reached at a time t* such that r'^'^et^ ~ 1 (for the argument 
in the exponential to be of order 1). Since tq ~ [tcou — ^Y^"^ in the scahng regime, we get 
tcoii — t* ~ etcoii- Therefore, the width of the scahng regime, 6t = t^oii — t*, behaves hke 

(85) 5t ~ Uoiie ~ e^'\ 

estabhshing v = 1/2. Close to the critical point, the collapse occurs at a very late time and the 
width of the scaling regime is very small. Therefore, if we are close to the critical point, it will 
be difficult to reach numerically the regime in which the results of sections |3.1| - |3.3| are valid. 

Regrouping all these results, and using again an analogy with critical phenomena, we expect 
that the central density obeys the following equation 

(86) p(0, t) = {t^u - t)-^G{t,oii{tcoii - t)), 

where tcoii ~ e~^^'^ and the scaling function G satisfies 

3 

(87) G'(O) = — , G{x) ~ pooVx, for x -^ +cx). 

2n 

4 Numerical simulations 

In this section, we perform direct numerical simulations of the SP system and compare the 
results of the simulations with the theoretical results of Sees. ^ and ^ In most of numerical 
experiments, we start from a homogeneous sphere with radius R and density p^, = SM/At^R^. 
This configuration has a potential energy Wq = —3GM^/5R. In the canonical ensemble the 
temperature is equal to T at any time. In the microcanonical ensemble, the initial temperature 
To is adjusted in order to have the desired value of A = 3/5 — 3RTq/2GM. By changing 
the temperature or the energy, we can explore the whole bifurcation diagram in parameter 
space and check the theoretical predictions of Sees. ^ and |[ In the numerical work, we use 
dimensionless variables so that M = i? = G = ^ = l. 

4.1 Microcanonical ensemble 

We first solve the SP system with the constraint (^ insuring the conservation of energy. We 
confirm the predictions of the thermodynamical approach in the microcanonical ensemble. For 
A = 0.334 < Ac, the quantities p(0,t), T(t), rxit) and S{t) converge to finite values and 
the system settles down to a stable thermodynamical equilibrium state with a density contrast 
TZ ~ 596 less than the critical value ~ 709 found by Antonov [0. At large distances, the density 
decays approximately as r~^ like the singular isothermal sphere [^. For A = 0.359 > A^ the 
behavior of the system is completely different: p(0,t) and T(t) diverge to +oo and rxit) goes 
to zero in a finite time tcoii- We were able to follow this "gravothermal catastrophe" up to a 
density contrast 7?. ~ 10^. The entropy S{t) also diverges to +cx3, but its evolution is slower 
(logarithmic). For A = 0.335 = A^, the system first tends to converge towards an equilibrium 
state but eventually collapses. 

In Fig. 1^, we plot the inverse of the central density as a function of time for different 
values of A. For short times, the density is approximately uniform, as it is initially. In that 
case, the diffusion term in Eq. (^ is negligible and the system evolves under the influence of 
the gravitational term alone. Using the Poisson equation (|^), the Smoluchowski equation (|^) 
reduces to 

, , dp AtxG n 

(««) i = — " ■ 
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Solving for p(t), we get 



^9) 



p(0,t) = p. 



^ AtiG 

1 + ^p* t + ... 



(t^O), 



where p^, is the initial density. Over longer time scales, a pressure gradient develops and the 
two terms in the right hand side of Eq. (^ must be taken into account. The system first reaches 
a plateau with density ~ poo (corresponding to an approximate balance between pressure and 
gravity) before gravitational collapse takes place eventually at t ~ tcoii- In Fig. |^, we see that 
the collapse time tcoii depends on the value of A and increases as we approach the critical value 
Ac. To be more quantitative, we plot in Fig. ^ the collapse time tcoii as a function of the distance 
to the critical point A — A^. A scaling law is observed with an exponent ~ —0.4 close to the 
predicted value —1/2 (see Sec. p^) . 
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Figure 2: Time evolution of the central density for different values of A. The central density 
p{0,t) becomes infinite in a finite time tcoii{A.) depending on the value of energy A (labeling 
the curves). The figure shows that the collapse time diverges as we approach the critical value 
Ac = 0.3345 for which a local entropy maximum exists. 

During the late stage of the collapse, the density profiles are self-similar that is, they differ 
only in normalization and scale (Fig. ^). Indeed, if we rescale the density by the central density 
and the radius by the King radius, the density profiles at various times fall on to the same 
curve (Fig. ^. The invariant profile is compared with the scaling profile f{x) corresponding 
to a = amax and the agreement is excellent, except in the tail. This small discrepancy can be 
ascribed to the next correction to scaling (see section |3]^) which generates a power law profile 
between ri and r2 with an index 7 < 2. We have checked that the logarithmic slope of the 
profile at r = i? is equal to —1] in agreement with the boundary condition (^). However, this 
relation only holds in a tiny portion of the curve (invisible in Fig. ^ so that the "effective 
slope" is more consistent with a value a ~ 2.2. In Fig. P, we plot the inverse central density 
as a function of time. It is seen that, for t -^ tcoii, the central density diverges with time 
like {tcoii — t)~^ in good agreement with the theoretical expectation. The slope of the curve 
in Fig. ^ is approximately —0.313 but is consistently getting closer to the theoretical value 
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Figure 3: Evolution of the collapse time tcoii with A. The figure displays a scaling law tcoii ~ 
(A — Ac)^^ with S ~ 0.4 close to the theoretical value 1/2. 

—1/5.178... ~ —0.193 corresponding to a = amax as A increases, or as t approaches tcoii (the 
small difference is attributed to non scaling corrections, as discussed in Sec. ^131) . Note that 
a value of a = 2 would yield a much larger slope — 27r/3 ~ —2.094 (see Sec. p.2|) , which is 
clearly not observed here. Therefore, the simulations are consistent with a value of a = amax, as 
expected on physical grounds. This value a = a^nax is also consistent with the slow but existing 
divergence of the temperature. Indeed, the slope of the curve in Fig. |^ is approximately —0.1 
in agreement with the theoretical expectation. 

To study the development of the instability for short times, we start from a point on the 
spiral of Fig. |l| close to Ac but with a density contrast TZ > 709 (we have taken A = 0.3344 
and TZ = 779). This isothermal sphere, with density profile Peq{r), is linearly unstable as it is a 
saddle point of entropy (see Sec. |2.3|) . In Fig. |^, we have represented the density perturbation 



profile 6p{r,t)/peq{r) = {p{r,t) — peq{r)) / peq{r) that develops for short times. This density 
profile presents a "core-halo" structure (i.e. it has two nodes) in excellent agreement with 
the stability analysis of Padmanabhan (we have computed the exact theoretical profile to 
compare quantitatively with the simulation). 



4.2 Canonical ensemble 

We now solve the SP system with a fixed temperature T. We confirm the results of the 
thermodynamic approach in the canonical ensemble. When t] < t]c the system converges to 
an equilibrium state while it collapses for i] > rjc (isothermal collapse). The collapse time tcoii 
scales with rj — rjc (see Fig. P) with an exponent ~ —0.6 close to the theoretical value —1/2. 

In Fig. ^, we plot the scaled density p{r,t)/p{0,t) as a function of the scaled distance 
r/rx{t) at different times. The curves tend to superimpose but the thickness of the line indicates 
that we do not have a strict self-similar regime (in agreement with our theoretical analysis). 



Indeed, the invariant profile f{x) computed in section |3.2| matches the numerics very well in 
the core but does not adequately describe the halo. The difference is due to the non scaling 
part F{r,t) that accounts for the mass conservation. In Figs. 111^, the result of the numerical 
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Figure 4: Evolution of the density profile for A = 0.359 > Ac at different times. Starting from 
a uniform distribution at t = 0, the system develops a "core-halo" structure with a shrinking 
core. From this figure, we may suspect that the evolution is self-similar, i.e. the density profiles 
at different times can be superimposed by an appropriate rescaling. 







-2 



/ — \ 

o 

Q. 


-4 








C 


-6 






Scaling profile 




-8 






(«=«max) 




-10 










6 


-4 


-2 










ln(r/g 



Figure 5: This figure represents the (quasi) invariant density profile obtained for A = 0.359 by 
rescaling the density by the central density and the radius by the King radius. It is compared 
with the theoretical profile f{x) calculated by solving Eq. (|33[) with a 
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Figure 6: Time evolution of the inverse central density for A = 0.359. This curve displays a 
scaling regime l/p(0, t) = Aitcoii — t). The slope of the curve A ^ —0.313 is of the same order 
as the theoretical value —1/5.178 = —0.193 corresponding to a = amax- The small deviation 
is attributed to non scaling corrections. 
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Figure 7: Time evolution of the temperature for A = 0.359. The curve displays a scalin; 
regime T ~ {tcoii — 1)~^ . The value of 7 ~ 0.1 is in agreement with the theoretical value (|6^ 

for a = amarr.. 
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Figure 8: First mode of instability in the microcanonical ensemble. The clean line is obtained 
by solving the eigenvalue equation (pSf ) with A = and the broken line is obtained from the 
numerical simulation of the SP system. The profile of density perturbation presents a "core- 
halo" structure. 
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Figure 9: Evolution of the collapse time tcoii with rj. The figure displays a scaling law tcou ~ 
iv — Vc)'" with u ~ 0.6 close to the theoretical value 1/2. 
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simulation is compared more precisely with the full theoretical prediction involving the non 
scaling term. The agreement is excellent throughout the whole domain. In the core, the profile 
is dominated by the scaling part which implies a r~^ behavior at moderately large distances. 
As explained previously and in Sec. ^.2| , this scaling behavior ceases to be valid near the wall 
and the contribution of the non scaling part is clearly visible. Its influence on the density 
profile remains weak but when the density is multiplied by r^, this non scaling profile has a non 
negligible contribution to the total mass. In Fig. ^, we see that the central density diverges 
with time as {tcoii — t)~^- The slope of the curve is approximately equal to 2 in good agreement 
with the theoretical prediction 2n/3 ~ 2.1 of section 
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Figure 10: Self-similar profile for rj = 2.75 > rjc- This (quasi) invariant profile is compared with 
the analytical scaling profile f{x) with a = 2. Deviation from the pure scaling law is due to 
non-scaling corrections that compensate for the excess of mass contained in the scaling profile. 



In Fig. [Ij, we study the early development of the instability for r] ~ ?7c. More specifically, 
we start the simulations from a point on the spiral of Fig. |1| with t] = 2.510 and 7^ = 42 > 32.1. 
This isothermal sphere is linearly unstable in the canonical ensemble as it is a saddle point 
of free energy (see Sec. |2.3| ) and the perturbation profile that develops for short times is 
shown in Fig. |l^. It is in excellent agreement with the first mode of instability calculated by 
Chavanis |]TD|] in the canonical ensemble. This profile does not present a "core-halo" structure, 
in contrast with the first mode of instability in the microcanonical situation. We have also 
plotted the perturbation profile for an isothermal sphere located near the second extremum 
of temperature [t] = 1.842...) at which a new mode of instability appears |]T0|. This second 
mode of instability has a core-halo structure (Fig. |15D. Of course, the perturbation profile that 



develops is a superposition of the first two modes of instability, but we see that its structure is 
dominated by the contribution of the second mode. 

In order to check the inequivalence of microcanonical and canonical ensembles in the region 
of negative specific heats, we started the simulation from an isothermal sphere with a density 
contrast comprised between 32.1 and 709. In the first experiment, the energy is kept fixed using 
the constraint (^. In that case, it is found that the sphere is linearly stable as it is a local entropy 
maximum. However, if the temperature is fixed instead of the energy, the sphere is now unstable 
as it is a saddle point of free energy. This clearly demonstrates in the framework of our simple 
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Figure 11: We plot the numerical finite-time density profile for rj = 2.75 (N.S.), at a time 
for which the central density is p(0,t) ~ 124.9 ~ 28.8poo- This is compared to the exact 
scaling profile Pof{r/ro) (Theory), with / given by Eq. (|38|) , and po = ^p(0,t) ~ 261.6 and 
ro = (r/po)-^/2 ^ 0.0373 (p(ro, t)/p(0, t) = 14/27 ^ 0.519). We also plot the asymptotic density 
profile, Pas = (vrr^r^)"^, valid for ro <^ r ^ 1. In this region, the correction to scaling is 
negligible. 



dynamical model that the microcanonical and canonical ensembles are not interchangeable for 
self-gravitating systems. This particular circumstance can be traced back to the non-extensivity 
of the system due to the long-range nature of the gravitational potential. This interesting 
problem is discussed in the review of Padmanabhan ||l|] and illustrated by Chavanis |^ for 
specific models of self-gravitating systems with a short-range cutoff (self-gravitating fermions 
and hard-spheres models). 

Since the stable isothermal configurations are only metastable (i.e., local maxima of a ther- 
modynamical potential), the value of energy or temperature is not sufficient to completely 
determine the evolution of the system: depending on the shape of the density profile, an initial 
configuration with A < Ac or r] < ?7c can either reach a quiescent equilibrium state or collapse. 
The actual evolution of the system depends whether the initial configuration lies in the "basin 
of attraction" of the local entropy maximum or not. Of course, the complete characterization 
of this basin of attraction is an impossibly complicated task because we would have to test 
all possible initial configurations. We have limited our study in the canonical ensemble to the 
case of unstable isothermal spheres located after the first turning point of temperature. These 
solutions correspond to saddle points of free energy. Therefore, a small perturbation (due here 
to numerical roundoff error) can destabilize the system and induce a dynamical evolution. The 
question is whether the system evolves towards the local maximum of free energy or undergoes 
gravitational collapse. Since we start from a saddle point of free energy, the two evolutions are 
possible depending on the form of the perturbation. In addition, depending on the location of 
the saddle point on the spiral (its density contrast), one of these evolutions may be preferred. 
The results of our study are displayed in Fig. |T6[ The isothermal spheres that experienced 
a complete collapse in our numerical experiments are marked with a symbol A while those 
that converged towards an equilibrium state are marked with a symbol •. A kind of structure 
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Figure 12: We plot the same numerical data (N.S.) as in Fig. |r^, but in the range 5ro < r < 1. 
This is compared with the theoretical density profile at t = tcoii obtained from Eqs. p8|) and 
(57). The parameters a ^ 5.0 and 6 ^ 5.1 are determined by maximizing p(l) (see text), 
although the full profile barely depends on a and b, as soon as b remains slightly greater than 
a, and b ~ 4.8 ~ 5.4. In this range, the theoretical profile is in excellent agreement with 
the numerical one. For instance, p(1)n.s. ~ 0.058 and p(l)Thcory ~ 0.057. In order to stress 
the quantitative agreement, we also plot the naive large r asymptotics of the scaling profile 
Pas = {iT'qr'^y^, for which p(l)as ~ 0.116. 




0.61 0.62 0.63 0.64 0.65 0.66 

t 

Figure 13: Time evolution of the inverse central density for rj = 3.5. This curve displays a 
scaling regime l/p(0,t) = B{tcoU — t). The slope S ~ 2 is close to the theoretical prediction 
27r/3 ~ 2.1. 
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Figure 14: First mode of instability in the canonical ensemble. The clean line is obtained 
by solving the eigenvalue equation ( p7\ ) with A = and the broken line is obtained from the 
numerical simulation of the SP system. The density profile does not present a "core-halo" 
structure. 
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Figure 15: Second mode of instability in the canonical ensemble. 
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Figure 16: Basin of attraction in the canonical ensemble. The isothermal spheres located after 
the first turning point of the spiral are unstable in the canonical ensemble. Depending on their 
position on the spiral (and the initial perturbation), they can either relax towards the local 
maximum of free energy with same temperature (•) or undergo a gravitational collapse (A). 

seems to emerge: it appears that the isothermal spheres undergoing gravitational collapse in 
the canonical ensemble are concentrated near the vertical tangent. We have found a similar 
structure in the microcanonical ensemble with a concentration of points undergoing gravita- 
tional collapse concentrated this time near the lower horizontal tangent. However, as indicated 
previously, this apparent structure is relevant at best in an average sense since other initial 
perturbations of the same saddle point may lead to a different evolution. In any case, these 
results confirm that the maxima of entropy or free energy are not global maxima since they 
do not attract all initial conditions. While homogeneous spheres with A < Ac and rj < rjc 
always seem to converge towards equilibrium, centrally concentrated systems with the same 
control parameters can develop a self-similar collapse leading to a finite time singularity. In 
fact, considering Fig. [l^ again, we see that the central concentration is not the only condition 
for collapse since there exists highly concentrated states that also converge towards the smooth 
equilibrium profile with low density contrast (in that case, the evolution corresponds to an 
"explosion"). Therefore, the basin of attraction of the metastable equilibrium states seems 
to have a highly non trivial structure. The nonlinear stability of a linearly stable isothermal 
sphere (located this time before the first turning point of energy or temperature) is also of 
interest. Since it is not a global entropy maximum it can be in principle destabilized by a 
finite amplitude perturbation. However, this perturbation is expected to be huge so that, in 
practice, the stability of the isothermal spheres with low density contrast is extremely robust. 



This suggests that these metastable states can be very long lived [g^, |2^, ^ and physically 
relevant in an astrophysical context. 
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5 Conclusion 

This paper has discussed the thermodynamics and the collapse of a system of self-gravitating 
Brownian particles in a high friction limit. This approximation considerably simplifies the 
problem since the evolution of the full distribution function /(r, v, t) is simply replaced by the 
evolution of its lowest moments. We showed that the Smoluchowski-Poisson system presents a 
rich variety of behaviors and displays interesting phase transitions between equilibrium states 
and collapsing states depending on the value of energy and temperature. When the two evo- 
lutions are possible, the choice depends on a complicated notion of basin of attraction. This 
simple model also illustrates dynamically the inequivalence of statistical ensembles for systems 
with long-range interactions. 

An extension of our study is to consider rotating systems with conservation of angular 



momentum. The SP system can be generalized to include rotation ||20[ and is interesting to 
study isothermal configurations that are not spherically symmetric. When spherical symmetry 
is broken, it is possible that the system will fragment in several clumps and that these clumps 
will themselves fragment in substructures. This may yield a hierarchy of structures fitting one 
into each other in a self-similar way as suggested by theoretical considerations pT| , |T^ . It would 



be of interest to investigate whether the SP system can display a process of fragmentation and 
exhibit a fractal behavior. Numerical simulations are under way. 

There exists a close analogy between the statistical mechanics of self-gravitating systems 
and two-dimensional vortices [^, ^, ^. Following the pioneering work of Onsager [^, there 



has been some attempts to describe vortices as maximum entropy structures, with possible 
applications to oceanic and atmospheric situations (e.g., Jupiter's Great Red Spot). The re- 
laxation towards the maximum entropy state is usually described by a Smoluchowski-Poisson 
system which analyzes the evolution of the vorticity in terms of a diffusion and a drift. The 
diffusion is due to the fluctuations of the velocity field and the drift to the inhomogeneity of 
the vorticity field |3^. The SP system can be deduced directly from the Liouville equation by 



using projection operator technics [^ or from a phenomenological maximum entropy produc- 



tion principle [^]. It is interesting to note that, for point vortices, the Fokker-Planck equation 
directly has the form of a Smoluchowski equation whereas for material particles this is true 
only in a high friction limit. This is because, for point vortices, the phase space coincides with 
the configuration space while for material particles it involves the positions and the velocities 
of the particles. 

The Smoluchowski-Poisson system also appears in the description of biological systems 
like bacterial populations [^. The diffusion is due to ordinary Brownian motion and the 



drift models a chemically directed movement (chemotactic flux) along a concentration gradient 
(of smell, infection, food,...). When the attractant concentration is itself proportional to the 
bacterial density, this results in a coupled system morphologically similar to the one studied in 
the present paper. The question that naturally emerges is whether this coupling can lead to 
an instability for bacterial populations similar to the gravitational collapse of self-gravitating 
systems. This possibility will be considered in a forthcoming paper in which we consider self- 
similar solutions of the Smoluchowski-Poisson equation for different systems in various space 



dimensions 26 . 
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A Analytical study of the scaling equation 

In this Appendix, we study analytically the scaling equation (|33|) . To that purpose, we rewrite 
it in an equivalent albeit more convenient form. Let us introduce the function 

(90) g{x) = 47r / f{x')x'^dx', 



in terms of which Eq. (|33D becomes 

(91) fix) + -fix) = \±{x'f'ix) + fix)gix)}. 

a x-^ ax 

Multiplying both sides of equation (0) by x^ and integrating the resulting expression between 
and X, we obtain 



(92) gix) = 4nx 



2 xfi^) - (^f'i.x) 
3 — a + Ana fix) 



From Eqs. ( pO]) and (|92|) , we can derive a nonlinear recursion relation satisfied by the coefficients 
an of the series expansion of fix) in powers of x^ (as / is an even function). Writing 



+ 00 
4:71 ^-^ 

n=0 



(93) /(-) = iE(-l 

we find 

_ 2n + a 1 ^^ CLpO-n-p 

^ ^ """+' ~ "2«(n + l)(2n + 3)''" ^ 2(n+l) ^ 2^+3 ' 

This recursion relation leads to the large n behavior of a„: 
(95) a„~8r(r2 + ^ )r" + «(0, 



where r is an unknown constant related to the inverse radius of convergence of the series. For 
a = 2, the asymptotics given by Eq. (p5| ) with r = 1/2 is an exact solution of the recursion 
relation (ISl), as can be checked by direct substitution. Using the identities 



1 i^e Orr 2 ^ 

(96) Pix) = . = y (-1) "rV", P'ix) = ^ = - y (-l)"nrV", 

71=0 n=0 

the series ( p3D can easily be resummed leading to Eq. (p8|). 

B The case of cold systems (T = 0) 

For T = 0, the core radius is not given by the King radius (^) which is zero by definition. 
We still assume however that po^o ~ ^i where a is unknown a priori. The equation for the 
invariant profile is then given by 

(97) fix) + -fix) = -Mf{x)9ix)), 

a x^ ax 
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where g{x) is defined by Eq. (pOl). Multiplying Eq. ([97| ) by Attx'^ and integrating from to x 
we obtain 

(98) (7(x) - '^ ' 



3 — a + ATTaf{x) 
Using the relation f{x) = g'{x)/A'Kx'^, the foregoing equation can be rewritten 

(99) (a — 3)g{x) + xg\x) = a^g\x)g{x). 

x^ 

Introducing the change of variables m = x^, we get 

(100) 3^ = i^^^. 

du u — ag 

A separation of the variables can be effected by the transformation g = uh, yielding 

(101) .]~''\ dh=--. 
^ ' h{3h-l) 3 u 

This equation is readily integrated leading to the implicit equation 

(102) g(x) = \i^- - g{x] 

where A is an integration constant. As g{x) is an odd analytical function, Eq. (|102|) first implies 
that g{x) ~ Y' ^° ^^^^ f^^^ ~ 'h' Combining with Eq. (|3^), this yields p(0, t) = -^{tcoii—t)~^. 
Then, inserting g{x) — ^ ~ 2;^ in Eq. ( p.02|) , we find that x^ ~ rf,5(i-a/3) ^ leading to a = 6/5. 
Note finally that the scaling profile defined by the implicit equation ( |102[ ) can be written in the 
parametric form 

(103) ^(^) = i^TT7' 9{x) = ls'/\ x = s'/'(^l + V'^' 

where the constant A has been incorporated in the expression of the core radius Tq. 

In fact, for T = 0, Eq. (^ can be solved analytically. Since the diffusion term vanishes, 
this equation describes a deterministic motion where the particles have a velocity u = — |V$ 



directly proportional to the gravitational force (see Sec. |2.1| ). This deterministic problem can 
be solved exactly by adapting the procedure followed by Penston |^ in his investigation of the 
collapse of cold self-gravitating gaseous spheres. Let us consider a particle located at r(0) = a 
at time t = 0. We denote by p(a) the average density inside the sphere of radius a. The total 
mass inside radius a can therefore be expressed as Ma = ^p{a)a^. At time t, this mass is now 
contained in the sphere of radius r = r(t), where r{t) is the position of the particle initially 
at r = a. Using the Gauss theorem, the motion of the particle is described by the first order 
differential equation 

dr IGMa 

(1°^) Tt = -5 V- 

This equation can be integrated explicitly to give 



(105) r = a{l-^p{a)t 
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Let us first discuss the case where the system is initially homogeneous with density p{a) = Pq. 
In that case, all the particles (whatever their initial position) arrive at r = at a time tcoii = 
$,/4:7iGpq defined as the collapse time for T = 0. This expression represents a lower bound 
(reached for rj — > +cx)) on the value of the collapse time tcouij]) studied in Sec. P^. During the 



evolution, the sphere remains homogeneous with radius, density and free energy evolving as 
(106) R{t) = R{1 - t/t,ouf'\ pit) = Po(l - t/t,our\ J{t) = ^^^(1 - t/t^oii)''^'. 



Note that the free energy diverges at t = tcou, unlike in Sec. ^]2[ These results can also be 
obtained directly from Eq. @ which reduces, for a uniform density, to 

dp ^fl ^^\ 1 ^^ 47rG 2 



(107) _^ = v^-pVtj=-pA* = — , 

where we have used the Poisson equation (|^) to get the last equality. 

We now suppose that, initially, p(a) has a smooth maximum at the center so that 

(108) p(a)=po(^l-^^, 

for sufficiently small a, where A is a constant. In that case, Eq. ( |105| ) giving the position at 
time t of the particle located at r = a at t = becomes 



1/3 

(109) 



i-t\' 



Ayt^oii, 

At t = tcou, the time at which the central density becomes infinite, it reduces to r = a^^'^ /A^^^. 
It is now straightforward to obtain the full density profile at t = tcoii- Since the mass contained 
between a and a + da at t = arrives between r and r + dr at time t, we have in full generality 

(110) p{a)ATTa^da = p{r, t)ATTr'^dr, 

or, for sufficiently small a, 

/-,-,-,-, , ^ _, .0? da _ a? da 

(111) py^^^) = PW—-^-Po—-^■ 

r ar r^ ar 

At t = t^oih we get 

(112) p{r,t,,u) = l-p,A'"'r-'l\ 

5 

We have therefore recovered that, for T = 0, the density profile decreases algebraically with an 
exponent a = 6/5. We now extend this analysis to a time r = tcoii—t just before the singularity 
arises. Considering the limit a —>■ and t —>■ 0, Eq. (|109|) can be expanded to lowest order as 

Then, Eq. ( |111| ) leads, after some reductions, to the density profile 
(114) p(r,t) = ^^. 
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The central density corresponds to r = 0, i.e. a = 0. According to Eq. (|114| ) it evolves with 
time as 

Pot coll _ ^ /^ _^_^-l 



(115) 
Therefoi 


'e, if we define 


p(0,t) = 


(116) 




5a tcoii 



T AnG 



{tcoll — t) . 



ro 



3^2 X 1/2 X ^ X5/6 



tcoll 



we can express the density profile in the parametric form 

P{r,t) _ 1 r _ ^1/2 A , 3 



^^^^^ p(o',t) 1 + .' ro(t) ' r + 5' 
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which is equivalent to Eq. ( |103|) . According to Eq. (|115|) and ( |116| ), we have the scaling laws 
^0 ~ (tcoll — tY^^i p(0)^o ~ l J^s^ before the singularity occurs. Setting F = p/p{0) and 
X = r/ro, we easily check that F{x) = 1 — x^ + ... for x — *> and F{x) ~ (3/5)^/^x^^/^ for 
X — > +00. This solves the problem for T = 0. Now, if the temperature T is very small but 
non-zero, we expect the present scaling to hold provided that tq ^ r^lt), where tq is defined in 
section |^. This leads to a cross-over core density pg above which the T 7^ scaling of section 
0| will prevail. The density p^ can be estimated by equating ro = (T/GpoY^'^ to tq ~ Po • 
The T 7^ scaling then prevails when the density becomes high enough, p^ ~ (T / G)^^^^ . 

C Connection between dynamical and thermodynamical 
stability 

Let p be a stationary solution of Eq. (^ and 6p a small perturbation around this solution. The 
first and second variations of temperature respecting the energy constraint (^) can be expressed 
as 

(118) ^M5T+ f 5p<^dh = 0, 



(119) -MS^T + - I 5p5^ dh = 0. 

The critical point p is a local entropy maximum provided that the second variations of entropy 

(120) ,^s_!M(^ + !^^^WM,3, 



4 r2 2 T 2j p 

are negative for any variations that conserve mass to first order. Let us now linearize Eq. (|^) 
around equilibrium and write the time dependence of the perturbation in the form 6p ~ e^*. 
We get 



(121) X6p = V 



^{STVp + TV6p + 6pV<^ + pV6<l>) 
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Multiplying both sides of Eq. ( |121| ) by 5p/ p, integrating by parts and using the equilibrium 
condition TV p + pV$ = 0, we obtain 

(122) A j ^-^ dh = - I -^{T\/5p + 5pV<^)i5TVp + TV5p + (5pV$ + pW5<^) £r. 

We now remark that the second order variations of the rate of entropy production (|TT]) are 
given by 

(123) 5'^S= I ^{8TVp + TV5p + 8pV^ + pV8^fd?v. 

J pTi 

We can therefore rewrite Eq. (|122|) in the form 

\j^^dh= -6'S + I ^{5TVp + pV5^) 

(124) x{5TVp + TV5p + 5pV^ + pV5^)dh. 

Using the equilibrium condition, the last term in Eq. ( [L24| ) is clearly the same as 

(125) - j ^{6TVp + TV6p + 5pV^ + pV(5$) ( ^V$ - ^V5$ j dh. 
Taking the time derivative of Eq. (^ and using Eq. (^ we have at each time 

(126) E = -Mf - I \{TVp + pV$)V<l' dh = 0. 

2 y 4 

The energy constraint ( |126| ) must be satisfied to first and second order. This yields: 



(127) f \{6TVp + TV6p + (5pV$ + pS/6^)V<^ dh = -M6f = -M\6T, 

(128) f jiSTVp + TV5p + 5pV^ + pV(5$) V(5$ dh = -M5^f = 3MX6'^T, 



where we have used Eqs. ( |118|) -( pJ^ ) to obtain the last equalities. Substituting these relations 
in Eq. (|2|), we get 



(129) A^,M!,3,,?M(|T!_3M^H-^^«. 



Comparing with Eq. (|120|) , we finally obtain 

(130) 6^S = 2X6^S. 

Since 6'^S > 0, see Eq. (|123D, the sign of A is the same as that of 6'^S. If p is a local entropy 



maximum, then 6'^S and consequently A are negative for any perturbation: the solution is 
linearly stable. Otherwise, we can find a perturbation for which 6'^S, and consequently A, 
are positive: the solution is linearly unstable. We can easily extend the relation ( |130|) to the 
canonical ensemble with J instead of S. We have found the same relation for other types of 
kinetic equations (Chavanis, in preparation), so its validity seems to be of a very wide scope. 
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